!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

module jm_mthd_mod;
use dbg_mod;
use HYP_2F1_MODULE;
implicit none;
private;

public test_jm_mthd_mod;

contains;!=====================================================================

subroutine test_jm_mthd_mod();
  character(LEN=64) :: name = "test_jm_mthd_mod";
  call dbg_in(name); 
  call test_hyp_2F1();
  call dbg_out(name); 
end subroutine

SUBROUTINE test_hyp_2F1();  
  real(DP) :: e, p, atomz, lambda;
  e = 0.5D0; !todo eng grid
  p = SQRT(2.D0 * e);
  atomz = 1.D0;
  lambda = 1.D0;
  call test_hyp_2F1_pz(p,  atomz, lambda, EPS14);   
  call test_hyp_2F1_pz(p, -atomz, lambda, EPS14);   

  atomz = 0.D0;
  call test_hyp_2F1_pz(p,  atomz, lambda, EPS12);   

  e = 1.D0; !todo eng grid
  p = SQRT(2.D0 * e);
  atomz = 1.D0;
  lambda = 1.1D0;
  call test_hyp_2F1_pz(p,  atomz, lambda, EPS12);   
  call test_hyp_2F1_pz(p, -atomz, lambda, EPS12);   
  
end SUBROUTINE

SUBROUTINE test_hyp_2F1_pz(p, atomz, lambda, maxErr);  
  real(DP), intent(IN) :: p, atomz, lambda;
  real(DP), intent(IN) :: maxErr;
  integer, parameter :: NMAX = 10, LMAX = 10;
  integer :: n, L;  
  complex(DP) :: a, b, c, z, f, HYP_2F1;
  real(DP) :: diff, t, dL, dn, TEST_2F1;
  real(DP) :: eta, x, theta;

  t = -atomz / p;
  eta = p / lambda;
  x = (eta*eta - 0.25D0) / (eta*eta + 0.25D0);
  theta = ACOS(x);
  do L = 0, LMAX
    dL = real(L, DP);
    call dbg_int("L=", L);   
    do n = 1, NMAX
      dn = real(n, DP);
      a = CMPLX(-dL, -t);
      b = CMPLX(dn+1.D0, 0.D0);
      c = CMPLX(dn + dL + 2.D0, -t);
      z = EXP(CMPLX(0.D0, -2.D0 * theta));
      call dbg_int("n=", n);   
      call dbg_cmplx("a=", a);   
      call dbg_cmplx("b=", b);   
      call dbg_cmplx("c=", c);   
      call dbg_cmplx("z=", z);   
      f = HYP_2F1(a, b, c, z);    
      call dbg_cmplx("f=", f);   
      diff = TEST_2F1(a, b, c, z, f);
      call dbg_dp("diff=", diff);   
      call dbg_assert_abs_diff_LE("err 2F1=", real(diff,DP), D0, maxErr);
    enddo
  end do
end SUBROUTINE

end module